% *** Read a fid file ***
% - Created Dat : 17 October, 2013
% - Modified Date : 17 October, 2013
% - Branch from : readfid(filename, channel)
% - Author : Sungjin Kim
% - Description : 
%     Reading Agilent fid file, and then rip off only the data block of each
%     channel, then saving the data to variables respectively.
%     The output has real and imaginary number both.
% - "c_args" options :
%     c_arg{1} = 1 : display header information in command line.
%     c_arg{1} = 2 : write header information onto disk.
%     c_arg{1} = 3 : nothing.
% - the position of header info : cell{1,1} file header
%                                 cell{1,2} file status
%                                 cell{2,#} block header
%                                 cell{3,#} block status
% *************************************************************************

function [cmplx_fid, headerInfo] = readfidn(filename, blockIdx, traceIdx, c_args)

%     s_fileHeader=32;  % the bytes of file header %
%     s_blockHeader=28; % the bytes of block header %
   
   %% Read fid file and save the binary information %%
    fileID = fopen(filename, 'r', 'b');
    
   %% The file header information %%
    % the file header information %
    nblocks   = fread(fileID,1,'int32');
    ntraces   = fread(fileID,1,'int32');
    np        = fread(fileID,1,'int32');
    ebytes    = fread(fileID,1,'int32');
    tbytes    = fread(fileID,1,'int32');
    bbytes    = fread(fileID,1,'int32');
    vers_id   = fread(fileID,1,'int16');
    status    = fread(fileID,1,'int16');
    nbheaders = fread(fileID,1,'int32');
    dataPoints = np/2;      % only real number data

    % Decode status information of file header%
    s_data    = bitget(status,1);
    s_spec    = bitget(status,2);
    s_32      = bitget(status,3);
    s_float   = bitget(status,4);
    s_complex = bitget(status,5);
    s_hypercomplex = bitget(status,6);
    s_acqpar = bitget(status,8);
    s_secnd = bitget(status,9);
    s_transf = bitget(status,10);
    s_np = bitget(status,11);
    s_nf = bitget(status,12);
    s_ni = bitget(status,13);
    s_ni2 = bitget(status,15); 
    
    % Header information of file and blocks %
    headerInfo{1, 1} = [nblocks, ntraces, np, ebytes, tbytes, bbytes, vers_id, status, nbheaders];  % file header
    headerInfo{1, 2} = [s_data, s_spec, s_32, s_float, s_complex, s_hypercomplex, s_acqpar, s_secnd, s_transf, ...
                                s_np, s_nf, s_ni, s_ni2]; % file header status
 
                            
    %% Header information of blocks %%
    f_rm = 0;
    for indx = 1:1:nblocks
        s_fseek = fseek(fileID, (indx-1)*bbytes+32, 'bof');
        % error check %
        if s_fseek == -1
            return
        end
        
        % the block header information %
        scale     = fread(fileID,1,'int16');
        bstatus   = fread(fileID,1,'int16');
        index     = fread(fileID,1,'int16');
        mode      = fread(fileID,1,'int16');
        ctcount   = fread(fileID,1,'int32');
        lpval     = fread(fileID,1,'float32');
        rpval     = fread(fileID,1,'float32');
        lvl       = fread(fileID,1,'float32');
        tlt       = fread(fileID,1,'float32');

        % Decode status information of block header%
        bs_data    = bitget(bstatus,1);
        bs_spec    = bitget(bstatus,2);
        bs_32      = bitget(bstatus,3);
        bs_float   = bitget(bstatus,4);
        bs_complex = bitget(bstatus,5);
        bs_hypercomplex = bitget(bstatus,6);
        bs_more_blocks = bitget(bstatus,8);
        bs_np_cmplx = bitget(bstatus,9);
        bs_nf_cmplx = bitget(bstatus,10);
        bs_ni_cmplx = bitget(bstatus,11);
        bs_ni2_cmplx = bitget(bstatus,12);
                
        headerInfo{2, indx} = [scale, bstatus, index, mode, ctcount, lpval, rpval, lvl, tlt];  % block header
        headerInfo{3, indx} = [bs_data, bs_spec, bs_32, bs_float, bs_complex, bs_hypercomplex, bs_more_blocks, ...
                                bs_np_cmplx, bs_nf_cmplx, bs_ni_cmplx, bs_ni2_cmplx];   % block header status
        
        % Display header infomation %
        if(c_args{1} > 0 && c_args{1} < 4)
                % create emty string array %
                str_Message='';

                % stack strings %
                if f_rm == 0
                    str_Message = [str_Message, sprintf('\n************** File Header Information **************\n')];
                    str_Message = [str_Message, sprintf('The number of blocks is %d\n', nblocks)];
                    str_Message = [str_Message, sprintf('The number of traces is %d\n', ntraces)];
                    str_Message = [str_Message, sprintf('The number of data points is %d\n', np)];
                    str_Message = [str_Message, sprintf('The bytes of each element is %d\n', ebytes)];
                    str_Message = [str_Message, sprintf('The bytes of each trace is %d\n', tbytes)];
                    str_Message = [str_Message, sprintf('The bytes of each block is %d\n', bbytes)];
                    str_Message = [str_Message, sprintf('The software version is %d\n', vers_id)];
                    str_Message = [str_Message, sprintf('The number of block headers per block is %d\n', nbheaders)];

                    str_Message = [str_Message, sprintf('\n  ************** Status Information in File Header **************\n')];
                    str_Message = [str_Message, sprintf('  The status of data (0 = no data, 1 = data) is %d\n', s_data)];
                    str_Message = [str_Message, sprintf('  The status of spectrum (0 = FID, 1 = spectrum) is %d\n', s_spec)];
                    str_Message = [str_Message, sprintf('  The status of 32 is %d\n    *if the status of float is equal to 0,\n     the status of 32 = 0 means 16-bit integer or\n     the status of 32 = 1, 32-bit integer.\n    *if the status is 1, all of things is ignored\n', s_32)];
                    str_Message = [str_Message, sprintf('  The status of float (0 = integer, 1 = float) is %d\n', s_float)];
                    str_Message = [str_Message, sprintf('  The status of complex (0 = real, 1 = complex) is %d\n', s_complex)];     
                    str_Message = [str_Message, sprintf('  The status of hypercomplex (= 1) is %d\n', s_hypercomplex)];
                    str_Message = [str_Message, sprintf('  The status of acquisition parmeters (= 1) is %d\n', s_acqpar)];
                    str_Message = [str_Message, sprintf('  The status of second FT (0 = first FT, 1 = second FT) is %d\n', s_secnd)];
                    str_Message = [str_Message, sprintf('  The status of transposed (reversed) (0 = reqular 1 = transposed) is %d\n', s_transf)];
                    str_Message = [str_Message, sprintf('  The status of np dimension activation (= 1) is %d\n', s_np)];
                    str_Message = [str_Message, sprintf('  The status of nf dimension activation (= 1) is %d\n', s_nf)];
                    str_Message = [str_Message, sprintf('  The status of ni dimension activation (= 1) is %d\n', s_ni)];
                    str_Message = [str_Message, sprintf('  The status of ni2 dimension activation (= 1) is %d\n', s_ni2)];
                    str_Message = [str_Message, sprintf('  ****************************************************************\n')];
                end

                str_Message = [str_Message, sprintf('\n************** %d Block Header Information **************\n', index)];
                str_Message = [str_Message, sprintf('The scaling factor is %d\n', scale)];
                str_Message = [str_Message, sprintf('The block index is %d\n', index)];  
                str_Message = [str_Message, sprintf('The mode of data in block is %d\n', mode)]; 
                str_Message = [str_Message, sprintf('The completed transient value(average) for FID is %d\n', ctcount)]; 
                str_Message = [str_Message, sprintf('The f2 (2D-f1) left phase in phasefile is %d\n', lpval)];
                str_Message = [str_Message, sprintf('The f2 (2D-f1) right phase in phasefile is %d\n', rpval)];
                str_Message = [str_Message, sprintf('The level drift correction is %d\n', lvl)];
                str_Message = [str_Message, sprintf('The tilt drift correction is %d\n', tlt)];

                str_Message = [str_Message, sprintf('  ************ Status Information in first block header ***********\n')];
                str_Message = [str_Message, sprintf('  The status of data (0 = no data, 1 = data) is %d\n', bs_data)];
                str_Message = [str_Message, sprintf('  The status of spectrum (0 = FID, 1 = spectrum) is %d\n', bs_spec)];
                str_Message = [str_Message, sprintf('  The status of 32 is %d\n    *if the status of float is equal to 0,\n     the status of 32 = 0 means 16-bit integer or\n     the status of 32 = 1, 32-bit integer.\n    *if the status is 1, all of things is ignored\n', bs_32)];
                str_Message = [str_Message, sprintf('  The status of float (0 = integer, 1 = float) is %d\n', bs_float)];
                str_Message = [str_Message, sprintf('  The status of complex (0 = real, 1 = complex) is %d\n', bs_complex)];
                str_Message = [str_Message, sprintf('  The status of hypercomplex (= 1) is %d\n', bs_hypercomplex)];
                str_Message = [str_Message, sprintf('  The status of more blocks (0 = absent 1 = present) is %d (*)\n', bs_more_blocks)];
                str_Message = [str_Message, sprintf('  The status of np complex (0 = real 1 = complex) is %d\n', bs_np_cmplx)];
                str_Message = [str_Message, sprintf('  The status of nf complex (0 = real 1 = complex) is %d\n', bs_nf_cmplx)];
                str_Message = [str_Message, sprintf('  The status of ni complex (0 = real 1 = complex) is %d\n', bs_ni_cmplx)];
                str_Message = [str_Message, sprintf('  The status of ni2 complex (0 = real 1 = complex) is %d\n', bs_ni2_cmplx)];
                str_Message = [str_Message, sprintf('  *****************************************************************\n')];

            % Display or write the string %

            if(c_args{1} == 1)
                % display string %
                fprintf('%s',str_Message);
            elseif (c_args{1} == 2)
                if f_rm == 0
                    delete('headerinfo.txt');
                end                    
                f_rm = f_rm + 1;
                % write a log on disk %
                fID = fopen('headerinfo.txt', 'a', 'n', 'UTF-8');
                fwrite(fID, str_Message);
                fclose(fID);
            end
        else
            fprintf('[ERROR] - c_arg{1} has a wrong value.\n');
            return;
        end % end if of display or writing %
    end % end loop of block header stacking %
    
    
   %% get complex data which are in a position %%
    % Make empty array to save data %
    realFidData=zeros(dataPoints, 1);
    imagFidData=zeros(dataPoints, 1);
    
    % Error checker %
    position = (blockIdx-1)*ntraces+traceIdx;
    if position > nblocks*ntraces 
        fprintf('[ERROR] - the excessive fid number (this data file has only %d fid) \n', nblocks*ntraces);
        return;
    elseif position < 1
        fprintf('[ERROR] - you cannot input a negative value.\n');
        return;
    end
        
    % Calculate a seeking position %
    seek_byteIdx = (blockIdx-1)*bbytes + (traceIdx-1)*tbytes + 32 + 28;  % remove file and block header
    
    % Seek the file to read data %
    s_fseek = fseek(fileID, seek_byteIdx, 'bof');
    
    % Check error of the fseek operation %
    if s_fseek == -1
        return;
    end    
       
    % reading data which is a pair consisted of real and complex number %
    for idx2 = 1 : dataPoints
        realFidData(idx2) = fread(fileID, 1, 'float32');
        imagFidData(idx2) = fread(fileID, 1, 'float32');    
    end
    
    cmplx_fid = complex(realFidData, imagFidData);
    
    %% closing file %%   
    fclose(fileID);    
    return;

end  % the last line of function %